// 最基础的武汉市水土流失

var geometry = ee.FeatureCollection("users/joeyqmf83/WH_shp");
Map.addLayer(geometry,{},"WH_shp")
//这一部分是处理云盖
function maskS2clouds(image) {
  var qa = image.select('QA60');
  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask).divide(10000);
}
var dataset_sentinel2=ee.ImageCollection("COPERNICUS/S2_SR")
                        .filterDate('2020-01-01','2020-12-31')
                        .filterBounds(geometry)//
                        .filterMetadata('CLOUDY_PIXEL_PERCENTAGE',"less_than",20)
                        .map(maskS2clouds)

var mid_image=dataset_sentinel2.median().clip(geometry);
var Temp = ee.Image()

//////////////////////////////// Normalization /////////////////////////////////////
function Normalization(image){
  var minMax = image.reduceRegion({
  reducer: ee.Reducer.minMax(),
  geometry: geometry,
  scale: 30,
  maxPixels: 10e9,
  tileScale: 16
}); 

var temp = ee.ImageCollection.fromImages(
  image.bandNames().map(function(name){
    name = ee.String(name);
    var band = image.select(name);
    return band.unitScale(ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))))
                // eventually multiply by 100 to get range 0-100
                //.multiply(100);
})).toBands().rename(image.bandNames()).clip(geometry);
return temp
}


///////////////////////////// Yellow ////////////////////////////////////
function Yellow(image)
{
  var Y=image.expression(
    "(B3+B4)/2",
    {
      "B4":image.select('B4'),
      "B3":image.select('B3')
    }
    ).rename('Yellow')
  return Y
}
///////////////////////////// NRI ////////////////////////////////////
function NRI(image)
{
  var N = image.expression(
    "B8A/B3",
    {
      "B8A":image.select('B8A'),
      "B3":image.select('B3')
    }
    ).rename('NRI')
  return N
}
///////////////////////////// NDSI ////////////////////////////////////
function NDSI(image)
{
  var N=image.expression(
    "(B11-B5)/(B11+B5)",
    {
      "B5":image.select('B5'),
      "B11":image.select('B11')
    }
    ).rename('NDSI')

  return N
}

///////////////////////////// NDVI ///////////////////////////////////
var NDVI=mid_image.normalizedDifference(["B8","B4"]).rename("NDVI")

///////////////////////////// Water Mask ///////////////////////////////////
var mask_water_NDVI=NDVI.updateMask(NDVI.gt(0))//mask water

///////////////////////////// FVC ///////////////////////////////////
function calFVC(BestVI,region,scale){
    var num = BestVI.reduceRegion({
      reducer:ee.Reducer.percentile([5,95]),
      geometry:region,
      scale:scale,
      maxPixels:1e13
    });
    var min = ee.Number(num.get("NDVI_p5"));
    var max = ee.Number(num.get("NDVI_p95"));
    //quantile and combine
    var greaterPart = BestVI.gt(max);
    var lessPart = BestVI.lt(min);
    var middlePart =ee.Image(1).subtract(greaterPart).subtract(lessPart);
    //calculate FVC
    var tempf1=BestVI.subtract(min).divide(max.subtract(min));
    var FVC=ee.Image(1).multiply(greaterPart).add(ee.Image(0).multiply(lessPart))
                      .add(tempf1.multiply(middlePart))
    return FVC.rename('FVC');
}
var FVC=calFVC(NDVI,geometry,10)


///////////////////////////// slope ///////////////////////////////////
var srtm = ee.Image("USGS/SRTMGL1_003");
var dem = ee.Image(srtm).clip(geometry);
var slope = ee.Terrain.slope(dem);

Temp = ee.Image(Yellow(mid_image))
Temp = Temp.addBands(NRI(mid_image))
Temp = Temp.addBands(NDSI(mid_image))
Temp = Temp.addBands(NDVI)
Temp = Temp.addBands(FVC)
Temp = Temp.addBands(slope)

Temp = Normalization(Temp)

////////////////////////////////计算水土流失指数 /////////////////////////////////////
function SEUFM(image)
{
  var S=image.expression(
    "(1-FVC)*(1-NRI)*(1-slope)*Yellow*NDSI",
    {
      "FVC":image.select('FVC'),
      "slope":image.select('slope'),
      "NDSI":image.select('NDSI'),
      "Yellow":image.select('Yellow'),
      "NRI":image.select('NRI')
    }
    ).rename('SEUFM')
  return S
}
Temp = Temp.addBands(SEUFM(Temp))
Temp = Normalization(Temp)

Map.addLayer(Temp.select('SEUFM').updateMask(NDVI.gt(0)), {min:0, max:0.6, palette:["green", "blue", "red"]}, "SEUFM");

// Map.addLayer(Temp.select('Yellow').updateMask(NDVI.gt(0)), {min:0, max:0.8, palette:["green", "blue", "red"]}, "YSI");

// Map.addLayer(Temp.select('NRI').updateMask(NDVI.gt(0)), {min:0, max:0.7, palette:["green", "blue", "red"]}, "NRI");

// Map.addLayer(Temp.select('NDSI').updateMask(NDVI.gt(0)), {min:0.2, max:0.85, palette:["green", "blue", "red"]}, "NDSI");

// Map.addLayer(Temp.select('FVC').updateMask(NDVI.gt(0)), {min:0, max:0.8, palette:["red", "blue", "green"]}, "FVC");

// Map.addLayer(Temp.select('slope').updateMask(NDVI.gt(0)), {min:0, max:0.8, palette:["green", "blue", "red"]}, "slope");

// Map.addLayer(Temp.select('NDVI').updateMask(NDVI.gt(0)), {min:0, max:0.8, palette:["red", "blue", "green"]}, "NDVI");